Unidad Práctica 5: Introducción a la Visualización de Redes¶

En este notebook exploraremos los viajes de la Encuesta Origen-Destino 2012 de Santiago utilizando visualizaciones de redes.

Google Colab¶

Para ejecutar este notebook en Colab, primero ejecuten la siguiente celda. Luego el notebook se reiniciará.

In [1]:

%load_ext autoreload
%autoreload 2

%load_ext autoreload %autoreload 2

In [2]:

try:
    import google.colab
    !pip uninstall matplotlib -y
    !pip install -q condacolab
    
    import condacolab
    condacolab.install_mambaforge()
except ModuleNotFoundError:
    pass

try: import google.colab !pip uninstall matplotlib -y !pip install -q condacolab import condacolab condacolab.install_mambaforge() except ModuleNotFoundError: pass

Una vez que se ha reiniciado puedes ejecutar el resto de las celdas.

In [3]:

try:
    import google.colab
    !git clone https://github.com/zorzalerrante/aves.git aves_git
    !mamba env update --name base --file aves_git/environment-colab.yml
except ModuleNotFoundError:
    pass

try: import google.colab !git clone https://github.com/zorzalerrante/aves.git aves_git !mamba env update --name base --file aves_git/environment-colab.yml except ModuleNotFoundError: pass

Preámbulo y Carga de Datos¶

In [4]:

import sys
from pathlib import Path

AVES_ROOT = Path("../..") if not "google.colab" in sys.modules else Path("aves_git")

EOD_PATH = AVES_ROOT / "data" / "external" / "EOD_STGO"
EOD_PATH

import sys from pathlib import Path AVES_ROOT = Path("../..") if not "google.colab" in sys.modules else Path("aves_git") EOD_PATH = AVES_ROOT / "data" / "external" / "EOD_STGO" EOD_PATH

Out[4]:

PosixPath('../../data/external/EOD_STGO')

In [5]:

import matplotlib as mpl
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import geopandas as gpd

from aves.data import eod, census
from aves.features.utils import normalize_rows

# esto configura la calidad de la imagen. dependerá de tu resolución. el valor por omisión es 80
mpl.rcParams["figure.dpi"] = 192
# esto depende de las fuentes que tengas instaladas en el sistema.
mpl.rcParams["font.family"] = "Fira Sans Extra Condensed"

import matplotlib as mpl import matplotlib.pyplot as plt import seaborn as sns import pandas as pd import geopandas as gpd from aves.data import eod, census from aves.features.utils import normalize_rows # esto configura la calidad de la imagen. dependerá de tu resolución. el valor por omisión es 80 mpl.rcParams["figure.dpi"] = 192 # esto depende de las fuentes que tengas instaladas en el sistema. mpl.rcParams["font.family"] = "Fira Sans Extra Condensed"

In [6]:

comunas = census.read_census_map('comuna', path=AVES_ROOT / "data" / "external" / "censo_2017_R13")

comunas = census.read_census_map('comuna', path=AVES_ROOT / "data" / "external" / "censo_2017_R13")

In [7]:

zones = gpd.read_file(AVES_ROOT / 'data' / 'processed' / 'scl_zonas_urbanas.json').set_index('ID')
zones.head()

zones = gpd.read_file(AVES_ROOT / 'data' / 'processed' / 'scl_zonas_urbanas.json').set_index('ID') zones.head()

Out[7]:

	AREA	Zona	Com	Comuna	REGION	NOM_REGION	PROVINCIA	NOM_PROVIN	COMUNA	NOM_COMUNA	URBANO	TIPO	NOM_CATEG	SHAPE_Leng	SHAPE_Area	area_m2	geometry
ID
103	837.7500	103.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	13105	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	4.575649e+05	POLYGON ((-70.65492 -33.55500, -70.65517 -33.5...
104	998.8125	104.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	13105	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.725462e+05	POLYGON ((-70.67599 -33.55542, -70.67452 -33.5...
106	286.2500	106.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	13105	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	2.693838e+06	POLYGON ((-70.67891 -33.55911, -70.68019 -33.5...
115	857.4375	115.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	13105	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.534193e+05	POLYGON ((-70.67625 -33.55513, -70.67653 -33.5...
116	853.9375	116.0	13105	El Bosque	13	REGIÓN METROPOLITANA DE SANTIAGO	131	SANTIAGO	13105	EL BOSQUE	EL BOSQUE	CAPITAL COMUNAL	CIUDAD	0.152123	0.001391	7.184305e+05	POLYGON ((-70.66564 -33.55260, -70.66309 -33.5...

In [8]:

viajes = eod.read_trips(EOD_PATH)

# descartamos sectores que no sean relevantes en los orígenes y destinos de los viajes
viajes = viajes[
    (viajes["SectorOrigen"] != "Exterior a RM")
    & (viajes["SectorDestino"] != "Exterior a RM")
    & (viajes["SectorOrigen"] != "Extensión Sur-Poniente")
    & (viajes["SectorDestino"] != "Extensión Sur-Poniente")
    & pd.notnull(viajes["SectorOrigen"])
    & pd.notnull(viajes["SectorDestino"])
]

print(len(viajes))

viajes = eod.read_trips(EOD_PATH) # descartamos sectores que no sean relevantes en los orígenes y destinos de los viajes viajes = viajes[ (viajes["SectorOrigen"] != "Exterior a RM") & (viajes["SectorDestino"] != "Exterior a RM") & (viajes["SectorOrigen"] != "Extensión Sur-Poniente") & (viajes["SectorDestino"] != "Extensión Sur-Poniente") & pd.notnull(viajes["SectorOrigen"]) & pd.notnull(viajes["SectorDestino"]) ] print(len(viajes))

89541

In [9]:

personas = eod.read_people(EOD_PATH)
viajes_persona = viajes.merge(personas)

personas = eod.read_people(EOD_PATH) viajes_persona = viajes.merge(personas)

In [10]:

viajes_persona["PesoLaboral"] = (
    viajes_persona["FactorLaboralNormal"] * viajes_persona["Factor_LaboralNormal"]
)

viajes_persona["PesoLaboral"] = ( viajes_persona["FactorLaboralNormal"] * viajes_persona["Factor_LaboralNormal"] )

¿Cómo se relacionan las comunas de acuerdo a los viajes entre ellas, por propósito?¶

Primero, preparemos el GeoDataFrame de comunas. Tenemos que quedarnos solo con las comunas que nos interesan, y tenemos que asegurarnos que tenga los mismos nombres que en el DataFrame de viajes.

In [11]:

#zones.head()

#zones.head()

Hacemos dos cosas:

• Como tenemos las zonas urbanas, filtramos el GeoDataFrame para quedarnos solamente con aquellas comunas que están en el DataFrame de zonas.
• Hacemos un diccionario de código de comuna -> nombre de comuna a partir de las zonas y lo aplicamos.

In [12]:

comunas_urbanas = comunas[comunas['COMUNA'].isin(zones['Com'].unique())].drop('NOM_COMUNA', axis=1).copy()
comunas_urbanas['NombreComuna'] = comunas_urbanas['COMUNA'].map(dict(zip(zones['Com'], zones['Comuna'])))
comunas_urbanas.plot(facecolor="none", edgecolor="#abacab")

comunas_urbanas = comunas[comunas['COMUNA'].isin(zones['Com'].unique())].drop('NOM_COMUNA', axis=1).copy() comunas_urbanas['NombreComuna'] = comunas_urbanas['COMUNA'].map(dict(zip(zones['Com'], zones['Comuna']))) comunas_urbanas.plot(facecolor="none", edgecolor="#abacab")

Out[12]:

<AxesSubplot:>

El mapa es demasiado grande, así que lo recortaremos utilizando las zonas que conocemos:

In [13]:

from aves.features.geo import clip_area_geodataframe

from aves.features.geo import clip_area_geodataframe

In [14]:

bounding_box = zones.total_bounds
bounding_box

bounding_box = zones.total_bounds bounding_box

Out[14]:

array([-70.995712  , -33.76479118, -70.44502298, -33.17216184])

In [15]:

comunas_urbanas = clip_area_geodataframe(comunas_urbanas, zones.total_bounds, buffer=0.05)

comunas_urbanas = clip_area_geodataframe(comunas_urbanas, zones.total_bounds, buffer=0.05)

/home/egraells/resources/aves/src/aves/features/geo.py:27: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  ).pipe(lambda x: x[x.geometry.area > 0])

Calculamos la lista de aristas de nuestra red, es decir, la cantidad de viajes de una comuna a otra. En este caso, lo haremos con los viajes al trabajo.

In [16]:

matriz = (
    viajes_persona[
        (viajes_persona["Proposito"] == "Al trabajo")
        & (viajes_persona["ComunaOrigen"].isin(comunas_urbanas["NombreComuna"]))
        & (viajes_persona["ComunaDestino"].isin(comunas_urbanas["NombreComuna"]))
    ]
    .groupby(["ComunaOrigen", "ComunaDestino"])
    .agg(n_viajes=("PesoLaboral", "sum"))
    .reset_index()
)

matriz.head()

matriz = ( viajes_persona[ (viajes_persona["Proposito"] == "Al trabajo") & (viajes_persona["ComunaOrigen"].isin(comunas_urbanas["NombreComuna"])) & (viajes_persona["ComunaDestino"].isin(comunas_urbanas["NombreComuna"])) ] .groupby(["ComunaOrigen", "ComunaDestino"]) .agg(n_viajes=("PesoLaboral", "sum")) .reset_index() ) matriz.head()

Out[16]:

	ComunaOrigen	ComunaDestino	n_viajes
0	Calera de Tango	Calera de Tango	2400.450674
1	Calera de Tango	Cerrillos	47.478444
2	Calera de Tango	Colina	0.000000
3	Calera de Tango	Estación Central	0.000000
4	Calera de Tango	La Cisterna	222.818975

Podemos convertir esta lista en una matriz de adyacencia. Veamos como luce esta matriz con el esquema adjacency_matrix. Como vimos en clase, utiliza la misma codificación visual que el heatmap de tablas, por lo que podemos usar seaborn.heatmap para visualizarla:

In [17]:

fig, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(
    matriz.set_index(["ComunaOrigen", "ComunaDestino"])["n_viajes"]
    .unstack(fill_value=0)
    .pipe(normalize_rows),
    cmap="inferno_r",
    linewidth=1,
)

fig, ax = plt.subplots(figsize=(12, 9)) sns.heatmap( matriz.set_index(["ComunaOrigen", "ComunaDestino"])["n_viajes"] .unstack(fill_value=0) .pipe(normalize_rows), cmap="inferno_r", linewidth=1, )

Out[17]:

<AxesSubplot:xlabel='ComunaDestino', ylabel='ComunaOrigen'>

También podemos utilizar el dataframe matriz como una lista de aristas que podemos visualizar con un gráfico NodeLink. A diferencia de la clase pasada, donde debíamos calcular la posición de cada nodo, al utilizar información geográfica los nodos ya tienen una posición.

In [18]:

from aves.models.network import Network

od_network = Network.from_edgelist(
    # graficamos los viajes más representativos
    matriz[matriz["n_viajes"] > matriz["n_viajes"].quantile(0.75)],
    source="ComunaOrigen",
    target="ComunaDestino",
    weight="n_viajes",
)

from aves.models.network import Network od_network = Network.from_edgelist( # graficamos los viajes más representativos matriz[matriz["n_viajes"] > matriz["n_viajes"].quantile(0.75)], source="ComunaOrigen", target="ComunaDestino", weight="n_viajes", )

In [19]:

matriz

matriz

Out[19]:

	ComunaOrigen	ComunaDestino	n_viajes
0	Calera de Tango	Calera de Tango	2400.450674
1	Calera de Tango	Cerrillos	47.478444
2	Calera de Tango	Colina	0.000000
3	Calera de Tango	Estación Central	0.000000
4	Calera de Tango	La Cisterna	222.818975
...	...	...	...
1097	Ñuñoa	San Miguel	586.746182
1098	Ñuñoa	San Ramón	908.799520
1099	Ñuñoa	Santiago	21003.071894
1100	Ñuñoa	Vitacura	1199.871414
1101	Ñuñoa	Ñuñoa	18793.426963

1102 rows × 3 columns

In [20]:

from aves.visualization.networks import NodeLink

nodelink = NodeLink(od_network)
nodelink.layout_nodes(method='geographical', geodataframe=comunas_urbanas, node_column='NombreComuna')

from aves.visualization.networks import NodeLink nodelink = NodeLink(od_network) nodelink.layout_nodes(method='geographical', geodataframe=comunas_urbanas, node_column='NombreComuna')

/home/egraells/resources/aves/src/aves/models/network/layouts.py:165: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  node_positions = positions_to_array(self.geodf.geometry.centroid)

En esta ocasión tenemos una red dirigida. Sabemos que los viajes van desde una comuna de origen hasta una comuna de destino. Por tanto, necesitamos una manera de identificar la dirección de las aristas. Una manera de hacerlo es pintarlas con un gradiente de color donde el azul representa el origen y el rojo representa el destino:

In [21]:

nodelink.set_edge_drawing('origin-destination')

nodelink.set_edge_drawing('origin-destination')

In [22]:

od_network.detect_communities(
    method="ranked", hierarchical_covariate_type="discrete-poisson"
)

od_network.detect_communities( method="ranked", hierarchical_covariate_type="discrete-poisson" )

In [23]:

nodelink.set_node_drawing("plain", weights='in_degree', categories='community')

nodelink.set_node_drawing("plain", weights='in_degree', categories='community')

In [24]:

from aves.visualization.figures import figure_from_geodataframe

fig, ax = figure_from_geodataframe(zones, height=7)

# contexto
zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0)
comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1)

nodelink.plot(ax, nodes=dict(palette='PuRd', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2)

ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)')

fig.tight_layout()

from aves.visualization.figures import figure_from_geodataframe fig, ax = figure_from_geodataframe(zones, height=7) # contexto zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0) comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1) nodelink.plot(ax, nodes=dict(palette='PuRd', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2) ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)') fig.tight_layout()

En comparación con la matriz de adyacencia, en esta visualización además de identificar las relaciones entre comunas podemos apreciar el contexto geográfico. Antes de elegir una de estas dos visualizaciones, debemos considerar lo siguiente:

• ¿Nos interesa conocer la relación geográfica entre orígenes y destinos? Por ej., ¿queremos saber si comunas vecinas se comportan similar?¿Nos interesa la distancia de los viajes? En este caso, node_link es una buena solución.
• ¿Necesitamos ver todas las aristas? Si es así, el gráfico de node_link podría ser inadecuado, ya que no podemos cambiar la posición de los nodos, ni podemos ver las aristas en las que el origen y destino son iguales. Debemos usar adjacency_matrix.

Ejercicio Propuesto: realizar el mismo análisis para distintos propósitos de viaje. Realizarlo a nivel de zonas (ver README de aves para un ejemplo con zonas).

In [25]:

matriz_zonas = (viajes_persona[(viajes_persona["Proposito"].isin(["Al trabajo", 'Al estudio'])) & (viajes_persona['ZonaOrigen'] != viajes_persona['ZonaDestino'])
                             & (viajes_persona['ZonaOrigen'].isin(zones.index))
                             & (viajes_persona['ZonaDestino'].isin(zones.index))]
                    .groupby(['ComunaOrigen', 'ZonaOrigen', 'ZonaDestino'])
                    .agg(n_viajes=('PesoLaboral', 'sum'))
                    .sort_values('n_viajes', ascending=False)
                    .assign(cumsum_viajes=lambda x: x['n_viajes'].cumsum())
                    .assign(cumsum_viajes=lambda x: x['cumsum_viajes'] / x['cumsum_viajes'].max())
                    .reset_index()
)

matriz_zonas.head()

matriz_zonas = (viajes_persona[(viajes_persona["Proposito"].isin(["Al trabajo", 'Al estudio'])) & (viajes_persona['ZonaOrigen'] != viajes_persona['ZonaDestino']) & (viajes_persona['ZonaOrigen'].isin(zones.index)) & (viajes_persona['ZonaDestino'].isin(zones.index))] .groupby(['ComunaOrigen', 'ZonaOrigen', 'ZonaDestino']) .agg(n_viajes=('PesoLaboral', 'sum')) .sort_values('n_viajes', ascending=False) .assign(cumsum_viajes=lambda x: x['n_viajes'].cumsum()) .assign(cumsum_viajes=lambda x: x['cumsum_viajes'] / x['cumsum_viajes'].max()) .reset_index() ) matriz_zonas.head()

Out[25]:

	ComunaOrigen	ZonaOrigen	ZonaDestino	n_viajes	cumsum_viajes
0	Providencia	497	328	6161.176000	0.001764
1	Colina	738	737	5582.101576	0.003361
2	Lampa	754	753	4797.111922	0.004734
3	San Bernardo	774	760	4509.053943	0.006025
4	Vitacura	682	666	4474.565492	0.007306

In [26]:

zone_od_network = Network.from_edgelist(
    matriz_zonas[matriz_zonas['cumsum_viajes'] <= 0.5], source="ZonaOrigen", target="ZonaDestino", weight="n_viajes"
)#.largest_connected_component(directed=True)
zone_od_network.network, zone_od_network.num_vertices, zone_od_network.num_edges

zone_od_network = Network.from_edgelist( matriz_zonas[matriz_zonas['cumsum_viajes'] <= 0.5], source="ZonaOrigen", target="ZonaDestino", weight="n_viajes" )#.largest_connected_component(directed=True) zone_od_network.network, zone_od_network.num_vertices, zone_od_network.num_edges

Out[26]:

(<Graph object, directed, with 687 vertices and 1918 edges, 1 internal vertex property, 1 internal edge property, at 0x7f1bc458d850>,
 687,
 1918)

In [27]:

merged_zones = zones.reset_index().dissolve('ID')

merged_zones = zones.reset_index().dissolve('ID')

In [28]:

zone_nodelink = NodeLink(zone_od_network)
zone_nodelink.layout_nodes(method="geographical", geodataframe=merged_zones)
zone_nodelink.set_node_drawing("plain", weights='in_degree')
zone_nodelink.set_edge_drawing(method="origin-destination")

zone_nodelink = NodeLink(zone_od_network) zone_nodelink.layout_nodes(method="geographical", geodataframe=merged_zones) zone_nodelink.set_node_drawing("plain", weights='in_degree') zone_nodelink.set_edge_drawing(method="origin-destination")

/home/egraells/resources/aves/src/aves/models/network/layouts.py:165: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  node_positions = positions_to_array(self.geodf.geometry.centroid)

In [29]:

fig, ax = figure_from_geodataframe(zones, height=7)

# contexto
zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0)
comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1)

zone_nodelink.plot(ax, nodes=dict(color='white', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2)

ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)')

fig.tight_layout()

fig, ax = figure_from_geodataframe(zones, height=7) # contexto zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0) comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1) zone_nodelink.plot(ax, nodes=dict(color='white', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2) ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)') fig.tight_layout()

In [30]:

zone_nodelink.bundle_edges(
     method="force-directed", K=1, S=0.005, I=6, compatibility_threshold=0.65, C=6
)

zone_nodelink.bundle_edges( method="force-directed", K=1, S=0.005, I=6, compatibility_threshold=0.65, C=6 )

Out[30]:

<aves.visualization.networks.fdeb.FDB at 0x7f1af4eb7e20>

In [31]:

fig, ax = figure_from_geodataframe(zones, height=9)

# contexto
zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0)
comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1)

zone_nodelink.plot(ax, nodes=dict(color='white', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2)

ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)')

fig.tight_layout()

fig, ax = figure_from_geodataframe(zones, height=9) # contexto zones.plot(ax=ax, facecolor='#efefef', edgecolor='white', zorder=0) comunas_urbanas.plot(ax=ax, facecolor='none', edgecolor='#abacab', zorder=1) zone_nodelink.plot(ax, nodes=dict(color='white', edgecolor='black', node_size=150, alpha=0.95), edges=dict(alpha=0.5), zorder=2) ax.set_title('Viajes al trabajo en Santiago (en días laborales, EOD 2012)') fig.tight_layout()